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Preceding  pip  Hank 


1  INTRODUCTION 

The  dynamic  analysis  of  framed  steel  structures  under  earthquake 
loading  would  appear  to  be  well  in  hand,  and  many  studies  and  appli¬ 
cable  computer  programs  are  readily  available  [1,2,3].*  In  terms  of 
building  configuration,  structural  elements,  and  applied  loading,  the 
programs  are  general  in  nature,  but  all  assume  that  material  behavior 
can  be  approximated  in  terms  of  moment-rotation  joint  relationships. 

A  simple  bilinear  or  Ramberg-Osgood  [2,4,5]  curve  fit  Is  usually 
adapted  to  model  this  material  resistance  function  and  appears  adequate 
for  most  engineering  applications.  Currently  accepted  design  criteria 
for  steel  structures  with  respect  to  lateral  displacements,  ductility 
ratios,  and  lateral  forces  have  evolved  from  studies  using  the  above 
computer  programs.  However,  the  complexity  of  the  load-displacement 
hysteresis  loops  for  frame  type  structures,  where  buckling,  plastic 
straining,  large  displacements  and  joint  slip  are  active,  can  be 
judged  from  the  survey  given  in  [6]  of  computations  and  experiments 
on  model  structures. 

The  current  study  presents  some  results  obtained  during  the 
course  of  a  continuing  research  project  on  the  magnitude  and  distri¬ 
bution  of  energy  dissipation  in  dynamically  loaded  structures.  A 
significant  point  of  this  research  is  that  the  behavior  of  the  struc¬ 
ture  is  found  by  means  of  a  material  model  with  properties  derived 
from  the  exact  hysteretic  stress-strain  behavior  of  A36  steel. 
Measurements  of  the  cyclic  behavior  of  A36  steel  were  carried  out  in 
[8],  and  the  analytical  model  was  developed  and  validated  in  [7]  for 
uniaxial  and  multiaxlal  stress  states.  By  means  of  this  model,  one 
can  incorporate  the  results  of  advanced  research  in  material  behavior 

*Bracketed  numbers  refer  to  reference  list. 


into  engineering  analysis.  In  all  cases  known  to  the  writers,  cyclic 
material  behavior  has  been  deduced  from  monotonic  test  curves  while, 
in  fact,  highly  accurate  results  on  cyclic  behavior  are  readily 
available  [9,10,11],  The  development  of  this  higher  order  material 
model  supplements  the  current  activity  in  the  cyclic  testing  of 
structures  through  large  strain  reversals  [12], 

Since  the  material  resistance  is  synthesized  from  knowledge  of  a 
polntwise  stress-strain  relation  rather  than  a  moment-rotation  rela¬ 
tion  at  a  joint,  the  finite  element  method  rather  than  a  stiffness 
method  is  used  to  model  the  structure.  The  element  used  is  a  simple 
beam-column  with  the  stiffness  matrix  obtained  by  numerical  integra¬ 
tion  along  the  length  and  through  the  depth  of  a  member  element.  A 
finite  element  plane  stress  analysis  of  a  simply  supported  beam, 
cycled  through  three  load  reversals  with  large  plastic  strains,  was 
used  in  [7]  to  check  the  material  model  under  discussion.  Further 
analysis  of  the  same  problem,  using  a  beam-column  element  and  ideal¬ 
ized  material  properties,  is  given  in  [13].  The  purpose  of  the  study 
in  [13]  was  an  assessment  of  the  accuracy  of  numerical  solution  tech¬ 
niques  with  a  particular  reference  to  cyclic  loading.  The  refined 
material  model  has  been  incorporated  in  a  general  purpose  finite 
element  computer  program  [14]  capable  of  analyzing  the  nonlinear 
static  and  dynamic  responses  of  engineering  structures. 

The  objective  of  this  research  is  to  formulate  an  analytical 
procedure  for  determining  overall  structural  energy  dissipation 
properties  on  the  basis  of  experimental  cyclic  material  behavior.  A 
projected  use  of  the  numerical  results  is  an  assessment  of  the  effect 
of  nonlinear  material  behavior  on  the  damping  properties  of  structures. 
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The  preliminary  results  given  here  are  aimed  at  verifying  the  overall 
analytical  model  and  the  accuracy  of  the  numerical  solution  procedures. 
For  the  case  of  a  simple  portal  frame  under  earthquake  excitation  it 
is  shown  that  considerable  differences  in  displacement  values  occur 
between  the  current  results  and  those  obtained  using  a  bilinear  kine¬ 
matic  hardening  stress-strain  relation. 

2  METHOD  OF  ANALYSIS 

(i)  Finite  Element  Formulation 

The  general  equation  of  motion  for  a  body,  derived  from  the  prin¬ 
ciple  of  virtual  work,  is  written  in  terms' of  initial  geometry  as 

[M]  {qt }  -  -  ;y  [B]T  {S}  dV  +  {Pt>  (1) 

where  [M]  is  a  consistent  or  diagonal  mass  matrix,  {q},  {q }  are  the 
generalized  displacement  and  acceleration  vectors.  {P}  is  a  vector  of 
equivalent  nodal  forces,  {s}  is  a  vector  of  generalized  stresses. 

The  matrix  [B]  transforms  generalized  displacement  increments  at  the 
nodes  to  generalized  strain  increments  at  any  point  in  an  element.  J.t 
is  defined  by  the  equation 

{AE}  -  [B]  {Aq>  (2) 

and  is  a  nonlinear  function  of  displacement.  The  expressions  used  in 
forming  [B]  in  the  current  study  are  given  in  the  Appendix. 

A  modified,  or  corrected,  linear  incremental  form  [14]  of  Eq.  (1) 
is  obtained  as 
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[M]  {Aqt+At}  -  -  $AB]T  {S>  dV  -  /v[B]T  {AS}  dV  +  {APt+At> 

+  (- [M]{q  }  -  /[B]T  {S}  dV  +  (P  })  (3) 

t  v  t 

The  incremental  stress  vector  {AS}  can  be  expressed  in  terms  of  {Aq} 
by  using  Eq.  (2)  and  the  following  linear  Incremental  stress-strain 
relation 


{AS}  *  [D  ]  {AE}  (4) 

ep 

The  construction  of  the  elastic-plastic  stress-strain  transformation 

matrix  [D  ]  with  reference  to  the  new  material  model  follows  the 
ep 

approach  outlined  in  [15] •  By  substituting  Eq.  (4)  in  Eq.  (3)  and 
following  procedures  described  In  [15],  Eq.  (3)  can  be  written  as 

[M|  WW  -  -  'V  (4W +  {fflW  +  {It>  (5) 

where  the  tangential  stiffness  matrix  [Kt]  is  formed  by  combining  the 
first  two  integrals  on  the  right  hand  side  of  Eq.  (3),  and  the  "load 
correction  term"  {lt)  is  the  contribution  of  the  terms  inside  the 
parenthesis  in  Eq.  (3).  {lt}  is  the  total  unbalanced  force  at  time  t. 
(ii)  Time  Integration 

The  controversy  over  an  optimum  temporal  integration  operator  has 
produced  a  large  volume  of  studies  in  recent  years.  Direct  methods 
[16]  are  available  for  estimating  the  accuracy  of  these  operators  when 
applied  to  linear  systems  [17],  but  it  seems  that  only  numerical 
experiments  can  provide  a  measure  of  their  worth  in  the  nonlinear 
case.  In  a  recent  analysis  [18]  of  the  performance  of  several  popular 
operators  (namely  the  Houbolt  [19],  the  Wilson  [20],  the  Newmark 
Beta  [21],  and  the  central  difference  operator),  it  was  noted  that  the 
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accuracy  and  stability  of  solutions  with  these  operators  were  very 
different  between  corresponding  linear  and  nonlinear  applications. 

The  overall  conclusion  arrived  at  in  [18]  Is  that  inaccuracies  in 
nonlinear  solutions  accumulate  very  rapidly  when  integration  operators 
are  used  with  time  inc*.  "ents  larger  than  that  of  the  stability  limit 
given  by  the  central  difference  operator. 

In  order  to  further  investigate  the  result  described  above,  both 
the  central  difference  operator  and  the  Newmark  operator  witli  8“l/4 
and  Y*l/2  were  used  in  the  current  study.  The  latter  operator  is 
generally  applied  in  earthquake  engineering  and  is  the  basis  of  the 
solution  method  in  [3].  The  Incremental  value  of  the  displacement 
vector  is  given  ,by  the  central  difference  operator  as 

{Aqt>  -  +  At;2  (6) 

The  solution  cycle  is  formed  by  calculating  the  displacement  increment 
from  Eq.  (6)  and  using  this  value  to  compute  the  integral  on  the  right 
hand  side  of  Eq.  (1).  The  acceleration  for  the  next  step  is  found  by 
solving  Eq.  (1)  for  {q‘t}.  The  operator  is  not  self-starting  and  suit¬ 
able  initial  values  are  outlined  in  [18], 

Because  of  its  implicit  form,  the  Newmark  operator  is  more  conven¬ 
iently  combined  with  a  linearized  equation  of  motion  such  as  given  by 
Eq.  (5).  The  increment  of  acceleration  is  given  by 

“VAt}  ■  ®VAtJ  -  4t  V  -  -T  'V’  (7) 

Substitution  of  Eq.  (7)  into  Eq.  (5)  yields 
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(^f1+<Kt»  {4W"  {SW  +  {It> 

+  (At{q  }  +  4*  f<ic»  (8) 

Ac  C  1  C 

In  this  approach  {Aqt+^}  is  first  calculated  from  Eq.  (8)  and  then 
substituted  back  Into  Eq.  (7)  to  find  the  acceleration  increment.  The 
velocity  vector  is  given  by 

<W  -  V  +  t  “V  ♦  <W>  m 

Damping  forces  of  the  equivalent  viscous  type  are  not  included  at 
present  since  they  would  complicate  the  Interpretation  of  the  effect 
of  material  nonlinearity  in  the  response. 

(ill)  Constitutive  Relation  and  Refined  Material  Model 
The  new  refined  material  model  is  an  extension  of  the  work  of 
Martin  [22]  and  Jhansale  [23]  and  is  described  in  detail  in  [7,8].  A 
mechanical  analog  of  the  basic  series  model  wit!  only  three  elements 
is  shown  in  Fig.  1.  E  ich  spring  of  the  system  has  a  linear  stress- 
strain  relationship.  The  springs  in  the  parallel  spring-slider  elements 
are  not  deformed  until  the  applied  stress  reaches  the  yield  stress  of 
the  slider  in  that  element.  If  the  applied  stress  in  an  element 
reaches  above  the  yield  stress,  then  the  difference  is  stored  in  the 
element  as  a  residual  stress.  The  tangent  modulus  E*  on  a  given  segment 
of  the  curve  is  obtained  from  the  relation 


where  the  summation  is  extended  over  all  those  spring-slider  elements 
which  have  yielded.  The  stress-strain  response  with  this  model  ..ollows 
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Figure  1.  Stress-strain  response  of  the  series  model. 
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Massing* s  hypothesis  where  the  closed  hysteresis  loops  are  the  same 
form  as  the  stabilized  initial  branch  OA  of  the  stress-strain 
curve  except  for  an  enlargement  by  a  factor  of  two  [10].  The  hypothe¬ 
sis  also  implies  that  after  a  small  hysteresis  loop  EFGH,  the  loading 
branch  will  follow  the  path  El  instead  of  EJ.  This  is  an  example  of 
the  "memory"  effect  which  is  created  by  the  distribution  of  residual 
stress  in  the  elastic  springs. 

This  basic  model  cannot,  however,  accurately  trace  the  hysteretic 
stress-strain  curve  of  A36  steel  for  two  reasons.  Firstly,  A36  steel 
in  its  "stable"  condition  does  not  obey  Massing's  hypothesis.  This  is 
shown  in  Fig.  2  where  the  cyclic  stress-strain  curve,  increased  by  a 
factor  of  two, does  not  accurately  describe  the  hysteresis  loop  shapes, 
and  the  upper  hysteresis  loop  tips  do  not  fall  along  the  upper  trace 
of  the  largest  loop.  Secondly,  A36  steel  cyclically  hardens  or 
softens  before  becoming  stable.  It  was  observed  In  [23]  and  [7]  that 
even  during  the  hardening  and  softening  process  each  upper  branch  of 
stress-strain  curve  (except  for  the  initial  flat  top  branch)  could  be 
fitted  to  the  double  skeleton  curve  for  the  "stable"  material  by  trans¬ 
lating  the  loop  along  the  elastic  line.  The  loop  is  translated  along 
the  elastic  line  until  the  upper  loop  tip  is  approximately  tangent  to 
the  double  skeleton  curve.  An  outer  trace  of  a  given  loop  is  then 
described  in  terms  of  this  skeleton  curve  and  an  appropriate  stress 

offset,  S  ,  which  defines  the  amount  that  the  loop  must  be  displaced 
os 

along  the  elastic  line.  In  Fig.  3  the  individual  branches  of  hystere¬ 
sis  loops  obtained  from  constant  amplitude  strain  controlled  tests 
have  been  fitted  to  the  double  skeleton  curve.  As  shown  in  the  figure, 
the  stress  offset,  Sog,  depends  on  the  number  of  reversals  and  the 


Figure  3.  Hysteresis  loop  branches  fitted  to  the  skeleton  curve 


17 


total  strain  range.  In  this  study  Sos  is  approximated  by  a  family  of 
linear  functions  of  the  number  of  reversals  and  the  total  strain 
range. 

The  analytical  material  model  also  incorporates  cycle-dependent 
mean  stress  as  described  in  [7],  but  this  feature  is  not  activated  in 
the  present  study.  An  added  feature  is  that  the  initial  response 
of  the  model  is  ietermined  by  a  set  of  constants  derived  from  the 
"flat-top"  stress-s .rain  curve  obtained  by  testing  a  virgin  specimen 
of  A36  steel.  The  spring  constants  and  slider  stress  values  for  sub¬ 
sequent  reversals  after  the  first  loading  are  found  separately  from 
the  skeleton  curve  obtained  from  constant  amplitude  strain  controlled 
cyclic  tests  of  steel  specimens.  The  procedure,  used  to  compute  the 
spring  constants  and  slider  stress  values, is  described  in  detail  in 
Appendix  A  of  [7].  Agreement  between  the  refined  material  model  and 
a  complicated  cyclic  stress-strain  history  for  a  virgin  steel  is  shown 
to  be  highly  satisfactory  by  the  results  depicted  in  Fig.  4. 

3  NUMERICAL  RESULTS 

(i)  Verification  of  Material  Model  in  Structural  Application 

In  order  to  test  the  validity  of  the  material  model  in  structural 
applications  the  beam  problem,  described  in  Fig.  5,  was  developed. 

Ten  spring-slider  elerf*nts  were  used  in  the  material  simulation,  and 
the  dimensions  of  the  beam  were  chosen  in  order  to  obtain  large 
plastic  strains  without  introducing  nonlinear  geometric  effects.  The 
beam  was  loaded  and  reloaded  through  three  cycles  as  shown  in  Fig.  3. 
In  the  earlier  investigation  of  Plummer  [7]  the  beam  half-span  was 
modeled  by  15  plane  stress  8-node  isoparametric  finite  elements  [24], 
With  a  view  to  later,  more  general  frame-type  structural  applications, 
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the  half-span  was  modeled  with  10  simple  beam-column  elements  In  the 
present  study.  The  generalized  strain-displacement  relations  on  which 
the  beam-coluion  element  is  based  are  outlined  in  the  Appendix. 

The  numerical  and  experimental  results  for  the  beam  are  compared 
in  Fig.  6.  The  geometric  shepes  of  both  curves  are  so  clearly  similar 
that  it  appears  that  a  slight  extension  of  the  elastic  range  would 
allow  perfect  agreement.  The  dimensions  of  the  beam  span,  cross- 
section,  and  concentrated  loading  also  make  it  likely  that  shear 
deformations  contribute  in  part  to  the  difference  in  results.  This  is 
investigated  in  Fig.  7  where  the  results  of  the  plane  stress  analysis 
from  [7]  are  compared  with  the  current  beam-column  results.  The 
continual  redistribution  of  the  two-dimensional  stress  state  Jn  the 
plane  stress  results  r^eans  that,  even  for  zero  hardening,  the  load 
will  continue  to  increase  after  the  critical  mid-section  has  become 
totally  plastic.  In  a  similar  example  Felippa  [25]  found  a  35  percent 
increase  in  load  capacity  beyond  the  limit  load  in  simple  bending  at  a 
displacement  of  approximately  five  times  the  elastic  limit  value.  The 
effect  in  this  case  is  to  give  the  plane  stress  curve  a  greater  slope 
on  the  initial  branch  which  leads  to  higher  residual  stresses  in  the 
spring  elements  before  the  first  reversal.  The  residual  stresses  for 
the  beam-column  case  are  practically  zero  except  for  one  element,  and 
this  means  a  flatter  response  will  occur  on  the  subsequent  reversal. 

It  is  felt  that  this  effect  rather  than  shear  displacement  causes 
most  of  the  difference  in  results  in  Fig.  7  at  the  end  of  the  second 
reversal. 

These  preliminary  results  indicate  that  further  tuning  of  the 
analytical  material  model  may  be  in  order.  Current  suggestions  in 
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this  respect  are  to  drop  the  initial  monotonic  "flat-top"  behavior  and 
use  the  cyclic  skeleton  curve  for  calculating  one  set  of  spring  and  slider 
values  for  all  cases.  Omission  of  the  "flat-ton"  first  cycle  may  result 
in  an  overall  improved  response.  It  was  suggested  in  Appendix  A  of  [7] 
that  fictitious  residual  stresses  be  built  in  at  the  end  of  the  "flat-top" 
cycle,  but  this  may  overly  complicate  the  parameters  of  the  model. 

The  constitutive  relations  for  the  plane  stress  analysis  in  [7]  are 
based  cn  Prager's  kinematic  hardening  theory.  This  theory  assumes  that  dur¬ 
ing  plastic  deformation  the  loading  surface  translates  as  a  rigid  body  in 
stress  space  maintaining  the  size,  shape  and  orientation  of  the  original 
yield  surface.  Plummer  [7]  satisfied  this  constraint  by  not  adding  the 
stress-offset  to  the  yield  stress  of  the  first  slider.  This,  however,  dis¬ 
torted  the  shape  of  the  double  skeleton  curve  represented  by  the  spring  slider 
model.  In  the  present  study  the  stress-offset  was  added  to  the  yield  stress 
of  the  first  slider  also.  A  general  combined  kinematic-isotropic  hardening 
theory  is  under  active  development  by  the  writers  with  a  view  to  relaxing  the 
restriction  of  constant  size  of  the  loading  surface. 

(ii)  Response  of  a  Portal  Frame  to  Seismic  Loading 
The  portal  frame  described  'n  [26],  was  modeled  with  three  finite 
elements  per  member  and  subjected  to  a  selected  four  seconds  (1.5  sec.- 
5.5  sec.)  of  the  El  Centro  NS  earthquake  acceleration  record  magnified 
by  a  factor  of  1.5.  The  equations  of  motion  must  be  adjusted  to 
reflect  this  type  of  loading  and  Eq.  (8)  is  replaced  by 

<i$1  +  I*t)>  (AV4t)  ’  {AW  +  fIt> 

2  2  (11) 

+  ^  <4‘  «t> +  T  <V  -  H  <V> 

where  {Aut+^t}  *s  now  the  generalized  displacement  relative  to  base  of 


rrf -v  ~i  tafert  i 


24 


the  structure.  The  vector  {Qt)  is  formed  by  the  components  of  base 

acceleration  input.  The  equilibrium  correction  term  is  also  changed  to 
ttt)  -  -  [M]({ut}  +  {Qt})  -  /V[B]T{S}  dV  +  (Pt)  (12) 

For  the  solution  with  the  central  difference  operator,  Eq.  (1)  is 
solved  as  it  is  written  for  the  absolute  acceleration,  and  the  rela¬ 
tive  structural  acceleration  is  obtained  from  the  relation 

(ut)  -  {qt)  -  {Qt>  (13) 

Displacements  relative  to  the  base,  for  the  purposes  of  calculating 
internal  forces,  are  found  in  the  usual  manner  as 

{Aut}  -  (A ut-At}  +  At2{ut_At}  (14) 

The  stiffness  matrix  [Kt ]  is  generated  numerically  using  three 
Gaussian  points  along  the  length,  and  the  cross-sections  of  the  I-beanis 
are  modeled  by  three  specially  weighted  integration  points,  one  placed 
in  each  flange  and  one  pt  the  center.  The  number  of  points  in  all 
cases  can  be  varied,  and  an  option  also  exists  which  allows  stresses 
to  be  computed  at  one  or  all  Gaussian  integration  points.  A  diagonal 
mass  matrix  is  formed  in  the  program  by  collapsing  all  rows  of  the 
consistently  formed  matrix  on  to  the  diagonal.  The  diagonal  mass 
matrix  is  used  with  Eq.  (1)  and  the  central  difference  operator,  and 
the  consistent  mass  matrix  is  used  in  Eq.  (11)  since  there  is  no 
apparent  computational  advantage  in  using  the  diagonal  form. 

The  dotted  line  in  Fig.  8  indicates  the  response  of  the  frame  with 
the  new  material  model  and  the  central  difference  operator.  In  the 
same  figure,  the  solid  line  passing  through  the  triangles  represents 


25 


the  response  with  kinematic  hardening  idealized  by  a  simple  bilinear 

6  3 

stress-strain  curve  given  by  E*=  30x10  psi,  ay*36xlOJ  psi,  and 

E’j 

—  ■  0.1.  A  comparison  of  these  two  curves  indicate  the  influence  of 
the  added  cyclic  hardening  effect  in  the  refined  material  model  over 
the  idealized  behavior.  The  results  clearly  demonstrate  a  cumulative 
effect  as  the  divergence  between  the  results  increases  with  each 
reversal  of  the  response.  It  is  emphasized  that  the  model  structure 
is  a  simple  portal  frame  loaded  with  only  a  fraction  of  a  typical 
earthquake  input,  and  greater  differences  can  be  expected  in  the 
analysis  of  more  complex  framed  structures  under  an  extended  earth¬ 
quake  input. 

The  solid  line  passing  through  the  circles  in  Fig.  8  shows  the 
values  given  by  DRAIN-2D  which  is  a  standard  stiffness-type  program 
developed  by  Kanaan  and  Powell  [3].  The  material  response  is  based  on 
a  bilinear  moment-rotation  joint,  or  node,  relationship,  and  the 
moment-axial  force  interaction  diagram  with  sh.:pe  code  2  [3]  was 
adopted.  Each  frame  member  was  divided  into  three  equal  elements,  and 
the  elastic  constants  and  hardening  modulus  mentioned  above  were 
employed  in  the  analysis.  Inspection  of  Fig.  8  shows  that  the  results 
with  DRAIN-2D  agree  in  general  with  those  found  using  bilinear  kine¬ 
matic  hardening.  The  divergence  of  these  curves  at  the  final  peak 
values  in  the  response  may  be  attributed  to  the  different  yield  cri¬ 
teria  and  to  the  manner  in  which  the  structural  stiffnesses  are  com¬ 
puted  in  the  respective  computer  programs.  The  agreement  between 
DRAIN-2D  and  the  results  of  the  refined  material  model  at  the  final 
peak  displacement  appears  to  be  fortuitous  since  the  displacement 
histories  are  considerably  different  up  to  this  point. 
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(111)  Selection  of  Time  Increment  for  the  Temporal  Integration 
Operators 

For  linear  problems  the  choice  of  an  implicit  versus  explicit 
integration  operator,  such  as  the  Newma^'.  over  the  central  difference 
operator,  is  often  governed  by  the  fact  that  the  latter  is  only 
conditionally  stable  with  respect  to  time  increment  size.  Also,  accu¬ 
rate  solutions  with  comparatively  large  time  increments  can  be  obtained 
from  Newmark's  operator  with  3sl/4  since  the  amplitude  of  each  mode  is 
conserved  [21].  However,  the  results  obtained  in  [18]  for  a  geometri¬ 
cally  nonlinear  simple  beam  problem  showed  that  all  implicit  methods 
suffered  serious  inaccuracies  for  time  increments  greater  than  that 
found  to  give  stable  results  with  the  central  difference  operator,  and 
that  the  Newmark  method  was  the  worst  in  this  respect. 

In  the  current  study  the  portal  frame,  described  in  Fig.  8,  was 
subjected  to  a  sinusoidal  base  excitation  in  order  to  find  suitable 
integration  time-steps  for  the  central  difference  and  the  Newmark 
operator.  The  response  wap  obtained  with  the  refined  material  model 
as  well  as  with  DRAIN-2D  which  also  uses  the  Newmark  operator,  and  the 
results  are  shown  in  Fig.  9.  The  limiting  stable  time  increment  size 
for  the  central  difference  operator  was  found  to  be  0.0025  seconds,  and 
an  apparently  stable  response  was  obtained  with  the  Newmark  operator 
for  a  time  step  of  0.00625.  Fig.  9  indicates  clearly  the  problem  of 
identifying  between  stable  but  not  necessarily  accurate  solutions.  In 
this  case  all  results  were  in  agreement  up  until  a  time  of  0.45  seconds 
when  both  DRAIN-2D  and  the  finite  element  solution  with  Newmark's 
method  begin  to  diverge  from  the  central  difference  solution.  There¬ 
fore,  it  was  decided  to  assess  this  phenomenon  more  accurately  by 
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running  the  problem  described  in  Fig.  8  with  decreasing  time  increments 
until  convergence  was  achieved  for  the  response  over  the  time  of  the 
earthquake  input.  The  results  obtained  with  DRAIN- 2D  for  Increments 
of  0.005  and  0.0025  are  shown  in  Fig.  10.  No  change  occurred  with  a 
reduction  in  step  size  to  .00125  and  so  convergence  Is  assumed  for  At 
equal  to  .0025.  The  problem  is  now  evident  in  that  apparent  convergence 
occurs  over  the  first  cycle  or  so  with  subsequent  slow,  but  finally 
appreciable,  divergence  after  approximately  two  seconds  of  the  response. 
The  recommendation  is  made  to  use  the  stable  central  difference  time 
increment  in  that  the  seemingly  apparent  stability  and  delayed  appear¬ 
ance  of  gross  inaccuracies  with  implicit  operators  could  lead  one  to 
believe  that  one  had  a  convergent  solution  at  larger  time  steps. 

The  results  under  discussion  in  this  section,  with  the  exception 
of  the  central  difference  case,  were  all  obtained  using  the  incremental 
form  of  the  equations  of  motion  with  an  elementary  equilibrium,  or  load, 
correction  term.  It  can  be  argued  that  an  iterative  solution  at  each 
step  would  improve  the  accuracy  of  the  implicit  solution  schemes.  This 
was  attempted  in  [18],  and  for  the  problem  studied,  the  iterated  scheme 
made  a  difference  only  at  those  time  increments  which  had  already  caused 
large  errors,  and  no  overall  significant  improvement  was  apparent. 
Further  study  in  this  respect  is  necessary,  but  the  economic  aspect 
of  an  explicit  central  difference  versus  an  iterated  implicit  solution 
may  be  the  deciding  factor  in  favor  of  the  former  method.  An  added 
consideration  is  the  simplicity  of  central  difference  solution  scheme 
and  its  considerable  advantage  with  respect  to  computer  storage  demand 


over  the  other  methods. 
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4  CONCLUSIONS 

An  analytical  procedure  has  been  described  which  incorporates 
constitutive  relations  based  on  a  higher  order  material  model.  A 
comparison  of  analytical  and  experimental  results  for  a  virgin  A36 
steel  coupon  and  for  a  simply  supported  beam  has  shown  the  capability 
of  the  refined  model  to  accurately  predict  the  cyclic  hysteretic 
material  behavior. 

The  responses  of  a  portal  frame  subjected  to  the  El  Centro  NS 
earthquake  acceleration  record  have  been  compared  for  the  refined 
material  model,  a  bilinear  kinematic  hardening  model,  and  the 
DRAIN- 2D  model.  The  significant  differences  in  the  computed  results 
have  clearly  demonstrated  the  effect  of  cyclic  hardening  on  seismic 
structural  response,  and  warrant  a  further  study  of  this  phenomenon 
with  respect  to  energy  dissipation  and  fatigue  life  of  steel  struc¬ 
tures.  The  preliminary  dynamic  results  suggest  that  for  comparable 
accuracy  the  time  increments  for  implicit  integration  operators  and 
for  the  central  difference  operator  must  be  of  the  same  order  in  the 
case  of  nonlinear  problems  with  complex  time-motion  histories. 
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APPENDIX 


BEAM-COLUMN  FINITE  ELEMENT 

The  axial  and  normal  displacements  U  and  V  at  a  point  (X,  Y)  of 
the  beam  are  approximated  by  a  linear  and  a  cubic  interpolation 
function  of  X  respectively. 


U 


ao 


+  a^X 


(A-l) 


V 


+  bjX  +  b2X2  +  b3X3 


(A- 2) 


where  a^,  ,  .  .  . ,  are  generalized  displacement  coefficients. 

It  is  convenient  to  separate  the  strain  E  at  (X,  Y)  into  its 
membrane  and  bending  components  E  and  E.  ,  where 


Eb  "  'Y  dx7 


From  Eqs.  (A-l)  and  (A-2) 


dU 

dX  "  al 


§  -  bl  +  2>>2X  +  3b3X2 


(A- 3) 


(A- 4) 


An  increment  of  strain  is  obtained  from  Eq.  (A— 3)  as 


dU  dU 
*dX  dX 


-A  ,  dV  /AdV. 
(idx>  +  dX  <4dX> 


(A-5) 
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Substituting  Eq.  (A-4)  into  (A-5) ,  AE^  and  AE^  can  be  expressed 
In  terms  c£  Increments  of  generalized  displacements  as 


'  >1 
AE 

m 

m 

0 

^x 

0 

dV 

dX 

»2£ 

^b 

0 

0 

0 

0 

-2Y 

-6XY 

l  J 

- 

- 
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> 


> 


(A-6) 


V  J 

From  Eq.  (A-6),  matrix  [B]  which  relates  generalized  Incremental 
strains  with  generalized  Incremental  displaceme  its  is  readily  identi- 


fled  as 

•* 

A 

1+du  o 

dV 

3X2  — 
dX 

U 

1+dX  ° 

dX 
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